The effect of thresholding on temporal avalanche 
statistics 



Lasse Laurson, Xavier Ilia, and Mikko J. Alava 

Department of Applied Physics, Helsinki University of Technology, FIN-02015 HUT, 
Finland 

E-mail: f irstname . lastname (at) tkk.fi 

Abstract. We discuss intermittent time series consisting of discrete bursts or 
avalanches separated by waiting or silent times. The short time correlations can 
be understood to follow from the properties of individual avalanches, while longer 
time correlations often present in such signals reflect correlations between triggerings 
of different avalanches. As one possible source of the latter kind of correlations in 
experimental time series, we consider the effect of a finite detection threshold, due to 
e.g. experimental noise that needs to be removed. To this end, we study a simple toy 
model of an avalanche, a random walk returning to the origin or a Brownian bridge, 
in the presence and absence of superimposed delta-correlated noise. We discuss the 
properties after thresholding of artificial timeseries obtained by mixing toy avalanches 
and waiting times from a Poisson process. Most of the resulting scalings for individual 
avalanches and the composite timeseries can be understood via random walk theory, 
except for the waiting time distributions when strong additional noise is added. Then, 
to compare with a more complicated case we study the Manna sandpile model of 
self-organized criticality, where some further complications appear. 
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1. Introduction 

Many systems in Nature are characterized by an intermittent avalanche- like response to 
slow externally applied driving pj]. Examples range from laboratory scale experiments 
on magnetic systems [2] and e.g. the sound emitted in a fracture test on a piece of 
paper [3J, to various applications in plasma physics and astrophysics [U El Ej, as well 
as in geophysics, including in particular earthquakes occurring due to the interaction of 
slowly moving tectonic plates [7j. Often the easiest way to characterize the dynamics 
of such systems is to record a global activity time series V(t), such as the acoustic 
emission (AE) amplitude in the paper fracture experiment [3], the induced voltage in a 
Barkhausen noise measurement [2] or the AE activity in martensites [8]. Such signals 
are often observed to be composed of apparently distinct bursts or pulses, which in the 
limit of slow driving are associated to distinct and typically spatially localized avalanches 
of activity occurring in the system. A typical feature of such bursts or avalanches is 
that the statistics of various measures associated to them appear to lack a characteristic 
scale, e.g. the avalanche sizes and often also durations are usually characterized by 
power law distributions. 

Another often made observation is that such avalanche- like signals appear to exhibit 
complex temporal correlations. On shorter time scales, the correlations can be thought 
to arise from individual avalanches, and a scaling theory relating the high frequency 
power spectra of such signals to the scaling properties of individual avalanches has been 
demonstrated to apply in a number of systems with avalanche dynamics [91 [Hll fill fT2]. 
Longer time correlations are often assumed to be due to correlations between the 
triggerings of different avalanches. Such apparently distinct bursts of activity are often 
found to be clustered in time. In the context of earthquakes, events are usually divided 
into large main events with smaller fore- and aftershocks occurring before and after the 
main shock, respectively. The rate of aftershocks is typically found to obey Omori's 
law, i.e. it decays as a power law in time after the main shock [13]. Similar conclusions 
have been obtained for the foreshocks [HI IT5] . 

Another possibility to detect correlations between different avalanches is to study 
the distributions of waiting times between two consecutive avalanches. If the avalanches 
are triggered by an uncorrelated process, one expects the waiting times to obey an 
exponential distribution due to the Poisson process -like random triggering of avalanches 
|16j . Consequently, deviations from this simple form can be interpreted to indicate 
the presence of correlations in the triggering process. In a number of systems, the 
distributions of waiting times between avalanches are found to be of power-law type 

PI HnH51 0312011211 122]. 

Such correlations can have various origins depending on the physical situation 
at hand, and it is often difficult to identify the mechanism operating in a specific 
experiment. In systems where the avalanches are triggered by the external driving 
(as opposed to situations like thermal creep in which thermal fluctuations trigger the 
avalanches), numerical studies show that the inter-avalanche correlations can be due 
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to the properties of the driving signal [HJ [23]. One must also be sure that a global 
activity signal is considered: In a local V(t)-signal, non-exponential waiting times can be 
observed due to the spatiotemporal fractal nature of critical avalanches [21] . In general, 
one can define a probability per unit time that an avalanche is triggered somewhere in 
the system. If other mechanisms for the observed correlations such as those mentioned 
above can be excluded, a possibility is that the occurrence of an avalanche somehow 
intrinsically affects this probability, thus leading to inter-avalanche correlations. 

In this paper we will focus on the effect of a finite detection threshold necessarily 
present in any real experimental situation. The presence of a background, either due to 
noise or processes coexisting with the intermittent avalanches in a measured activity time 
series V(t), forces one to apply a finite threshold level V t h, and define avalanches as the 
bursts exceeding this threshold [H] [12] . While such ideas have attracted some attention 
in the literature [25j EH EI], neither theoretically nor experimentally general studies 
have been carried out on the effect of thresholding on the avalanche statistics. Since 
an avalanche can be defined as a correlated sequence of activity, breaking the avalanche 
into smaller parts by thresholding, the ensuing "subavalanches" will be temporally 
correlated, being part of the same underlying avalanche. As an attempt to clarify 
the issues associated to this effect of thresholding, we consider a simple toy model of 
the signal V(t) corresponding to a single avalanche, namely an excursion of a random 
walk from the origin. Such a simple model is sufficient to demonstrate how e.g. the 
distribution of waiting times arising from the thresholding process assumes a power law 
form. We also discuss the same phenomenology in the stochastic Manna sandpile model 
of self-organized criticality [29]. We also consider the addition of external noise, and 
the properties of composite signals made of toy avalanches and waiting times from a 
distribution chosen a priori: how thresholding affects them with or without the noise. 
For avalanche sizes and power spectra, there are no real complications, but for the 
waiting times we find a mixture of effects including an apparent power-law regime with 
a different exponent, when noise is added. The paper is organized as follows: In the next 
Section, excursions of random walks as a model of an avalanche signal is considered. 
The Manna sandpile model is briefly discussed in Section 3. Section 4 finishes the paper 
with conclusions and a summary of future prospects. 

2. The random walk model 

As a conveniently simple toy model illustrating the mechanisms associated to the 
thresholding process of avalanche time series, we consider an excursion of a discrete 
random walk, obeying 



where rj is white noise with a bimodal distribution P(rj) = l/2(d v i + b~ v ,-i) an d time is 
updated in steps of magnitude At = 1. The excursion starts at t — from x = 0, and 
returns to the origin for the first time at t = T. Such a random walk bridge is then 
taken to model the V(t)-signal corresponding to a single avalanche of duration T. 



d t x = T], 
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Figure 1. An example of subavalanches obtained when a random walk bridge of given 
duration T is thresholded with a threshold level Vth , with the definitions of the different 
waiting times (or quiet times) Tq, fi and ff. 



The statistical properties of such excursions are well known pTTj 132"! 133] . The first 
return times, or avalanche durations T obey a power law distribution 

P(T) ~ T~ TT (2) 

for T ^> 1, with Ty = 3/2. Similarly, the distribution of avalanche sizes s = Jq V(t)dt 
is given by 

P(s) ~ (3) 
with r s = 4/3. The average shape of the excursion is given by a semicircle, 

V(t,T)=T^- 1 f shape (t/T), (4) 



where f s ha P e(x) = J8/7rJx(l — x) and 7 si = 3/2. The exponent j s t relates the average 



'7st 



avalanche size, ( S (T)) = (/ J V(t)dt) to the duration T as (s(T)) ~ T 
2.1. Random walk bridges of a given duration 

First we will focus on the statistical properties of the avalanches which emerge when 
random walk bridges of fixed duration T are thresholded as is showed if Fig. HJ 

The first quantity of interest is the number of avalanches N(T, Vth) as a function 
of the random walk bridge duration T and the threshold level V t h- The upper panel of 
Fig. [2] shows the number avalanches observed as a function of the threshold level Vth, 
for various durations T of the unthresholded excursions. These can be collapsed onto a 
single curve by using the ansatz 

N(T,V th )=T a f N (V th /T 1 / 2 ). (5) 

The value of the exponent a = 0.51 ± 0.01 is observed to be close to 1/2, a result that 
can be understood to follow from the known scaling of the number of zero crossings 
N zc (t) of a random walk function of time, N zc (t) ~ t 1/2 [33J. 

The lower panel of Fig. [2] displays the distributions of the maxima of the random 
walk excursions. Again, a good data collapse is obtained by using the scaling 

P(x max ,T) = T~^f Xmax (x max /T 1 ^), (6) 
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Figure 2. Top: The rescaled number of waiting times as a function of the threshold 
level Vth, for different avalanche durations T. Bottom: The rescaled distributions of 
the maximum values of the avalanche signal as a function of threshold level Vth, for 
different avalanche durations T. The dashed line is a fit to a log-normal distribution 
f(x; fJ,, a) = [1/(x(jV2tt)} exp(-(log(x) - (i) 2 / (2a 2 )) , with fx = 0.21 and a = 0.22. 

with an apparently universal scaling function f Xmax (y) that can be fitted well with a log- 
normal distribution. The two first moments of this fitted distribution are in agreement 
with the ones computed in Ref. [28] using the exact expression. 

When a non-zero threshold Vth is applied, one can define three different kinds 
of waiting times: the initial waiting time To, the intra-avalanche waiting times fj, 
and the final waiting time ff, see Fig. [TJ Notice that in this paper we denote all 
quantities x defined with a non-zero threshold level V t h like x, while x stands for 
the same thing with Vth = 0. Since for models with a symmetrical avalanche shape 
P(ff] T, V t h) = P(to] T, V t h), we focus on the distributions of f and fj. The former one 
can be collapsed using the scaling ansatz 

P(f ;T,Vth)~V th 2 g(f /V t l). (7) 

Due to the Markovian and translation invariance properties of the random walk, 
P(fi] T, Vth) is expected to be of a power law form with the same exponent as P(T) 
just the durations of the excursions below the threshold Vth), i- e - 

p(n; t, Vh) ~ f^fin/vDhin/T). (8) 

The cut-off scaling follows from the restriction that such excursions below the threshold 
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Figure 3. Scaling functions for the initial waiting time distribution (top) and for the 
intra-avalanche waiting time distribution (bottom). Different colours (online version) 
account for different values of r = Vth/T 1 / 2 = 0.3,0.4,0.5,0.6. For larger r-values the 
intra-avalanche waiting times fj do not obey Eq. ([9]) as for very large r the cut-off due 
to the finite avalanche duration becomes important, see Eq. (|8j). 



are bounded between and V t h, so that a maximum waiting time scaling as Vf h ensues, 
see Eq. fl4]). For random walk bridges of a fixed duration T, another constraint is that 
fi < T. However, in practice the cut-off scale due to the finite threshold value is reached 
first, and the scaling function fr(x) can be regarded as a constant. Thus, one can write 
a T-independent form for the distribution of internal waiting times, 

Pin; T, Vth) ~ P(n; V th ) ~ n' i/2 f{n/V? h ). (9) 

Fig. [3] displays data collapses of distributions of tq and fj. The intra-avalanche waiting 
times fj obey Eq. (Q only for r < 0.7 as for very large r the cut-off due to the finite 
avalanche duration becomes important, see Eq. (jSJ). 

Next we consider both the initial and the internal waiting times together, by defining 
f = fo U fj. The total distribution of waiting times f is then given by 

= P(f ;r,^) + (7V(r,V, t )-l)P(r i; r,V, t ) 

N(T,V th ) V ; 

For T>1, the number of subavalanches N(T, V t h) 3> 1, see Eq. (|3j) and Fig. [2j In this 
limit N(T, V t h) — 1 — N(T, V t h) and one obtains the total distribution of waiting times 
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Figure 4. Scaling functions for the total waiting time distribution for different values 
of r. 



for a fixed duration T and threshold level Vth'- 
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Now, using the more natural variables r 
equation as 

g(q/r 2 



T- a + f~^f(f/V t l). (11) 
one can rewrite the above 



and q 



T ' 



P(g; T, r) = T~ 3/2 



+ g- 3/2 / (q/r 2 ) 



(12) 



7 2 /Af( r ) 

which allows us to collapse the total waiting time distributions for a given value of r 
for all durations T, as can be seen in Fig. HI Moreover, the use of these variables also 
allows us to write the mean value of the total waiting time distribution as 

,2 



(f) = T 1 ' 2 



Mr) 



+ I 2 r 



(13) 



where h = fd£ &(0 and I 3 = fd£ r 1/2 /(0- 



Entire signal 

In order to generate an artificial signal, one can insert exponentially distributed waiting 
times (corresponding to uncorrelated triggerings of the avalanches) between the original 
(non-thresholded) avalanches and study the evolution of the waiting time distribution 
as the threshold level Vth is increased from zero as can be seen in Fig. [5] (top). Also the 
effect of adding white noise to the random walk signal is considered as can be seen in 
Fig. [5] (bottom). 

The total statistics arising from the intra-avalanche properties of an ensemble of 
excursions with duration distribution given by Eq. (j2J) at a fixed threshold is obtained 
as a convolution over P(T). For instance, one can write for the distribution of internal 
waiting times 

P{n; = J P{T)P{n; T, V th )dT ~ f^'fin/Vtl), (14) 
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Figure 5. Top: An example of a signal V^(£) consisting of two random walk excursions 
from the origin, showing the definitions of the avalanche duration T and the waiting 
times (or quiet times) f after thresholding the signal with a threshold level Vth indicated 
in the figure by the horizontal line. Bottom: The same signal as above, but with 
superimposed Gaussian white noise. The presence of noise forces one to apply a finite 
threshold, but affects also the thresholded avalanches. 



where the last step follows from the fact that P{fi, Vth,T) is independent of T. Similar 
considerations apply for the duration and size distributions of the subavalanches induced 
by the thresholding process, such that the thresholded avalanche durations T obey 
P(T) ~ T~ 3 / 2 . Similarly, thresholded avalanche sizes s 

= fo[V(t) - V th ]dt follow the 

scaling P(s) ~ s~ 4//3 , see Fig. [6j 

When considering the case with exponentially distributed waiting times r inserted 
between subsequent random walk excursions, and studying the distributions of all 
waiting times f = fj U f e (see Fig. EJ) as a function of the applied threshold Vth, an 
evolution from an exponential to a power law distribution is observed as the threshold 
V t h is increased, see Fig. [7J 

From an experimental point of view, all V r (t)-signals are noisy. In order to treat 
these real signals, one should add a minimum threshold to get rid of the background 
noise but keeping the maximum amount of information. Therefore, to study this effect 
in our artificial signals, we add a Gaussian white noise £ with mean zero and standard 
deviation a to the ^(t)-signal and consider different thresholds ranging from Vth = 4a 
to V t h = 8cr. Surprisingly, the exponent appears to change from t w = 3/2 to r w w 2, 
see Fig. |HJ This change is because the original (without noise) waiting times are broken 
into shorter ones by the noise, but we have not managed yet to find a mathematical 
derivation of such an exponent. 

These results demonstrate how a symmetry between the avalanche durations and 
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Figure 7. The probability distributions of waiting times (or quiet times) f = t% U f e 
(see Fig. [5]) between avalanches, for different threshold values, and waiting times 
between the non-thresholded excursions drawn from a probability distribution P (t) = 
(1/6) exp(— r/b), with b — 100. For Vth — one recovers the exponential P (t) 
distribution, while for large enough threshold values the distributions become power 
laws with an exponent t w « 1.5. 
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Figure 8. The probability distributions of waiting times (or quiet times) f between 
avalanches with additive Gaussian white noise, for different threshold values Vth and 
noise strengths a . The waiting times between the non-thresholded excursions drawn 
from a probability distribution Po(t). The power law exponent t w changes from 
the noise free value of t w = 3/2 to t w sa 2.0. Considering the same signal without 
the artificial waiting times (not shown) produces similar distributions, with the only 
difference being a decrease in the magnitude of the "bump" in the cut-off due to the 
artificial waiting times. 




Um 1 Vth 



x 10 



Figure 9. An example of a signal V(t) consisting of avalanches from the 2d Manna 
model with L — 1024. The horizontal line corresponds to a threshold level with 
V th = 160. 



quiet time intervals could give a natural explanation to observations of power law 
distributed waiting times in various experimental situations in which any kind of 
thresholding process is applied. To see how such considerations can be generalized 
to more realistic systems with non-trivial dynamics, we consider in the following the 
same phenomenology in the two-dimensional stochastic Manna sandpile model of self- 
organized criticality. 

3. Manna model 

Out of many possible model systems, we consider here the Manna sandpile model 
of self-organized criticality [29]. This stochastic model has the advantage over some 
other models that it is known to obey simple scaling, whereas models such as the 
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original sandpile model introduced by Bak, Tang and Wiesenfeld (BTW) J3D] has been 
observed to exhibit multiscaling [31]. In Ref. [26], the effect of thresholding on the 
statistics of quiet time intervals and avalanche durations was considered. Similarly to our 
observations above concerning the simple random walk model, the exponents describing 
the distributions of the avalanche durations and quiet time intervals were observed be 
the same, but with a value close to 5/3, different from both the random walk value of 3/2 
and the value found for the avalanche durations in the BTW model without thresholding 
|35j. One should notice, however, that the latter observations regarding the exponent 
values can be masked by the multiscaling exhibited by the BTW model [34|. 

The stochastic Manna model is defined on a ci-dimensional hypercubic lattice of a 
linear size L, with an integer variable zi assigned to each lattice site i. A site topples if 
the local variable reaches or exceeds a critical value z c = 2. In the toppling process, two 
grains will be redistributed from the toppling site to its two randomly chosen nearest 
neighbours. During a single time step, all the sites will be checked and those with 
Zi > z c = 2 will be toppled in parallel. Such a toppling can then cause one of the 
neighbouring sites of the toppling site to topple during the next time step. Thus, 
an avalanche of activity can be triggered from a single initial toppling event. Such 
triggerings happen due to the addition of new grains to random locations in the system 
at a slow rate. This external driving is balanced by allowing grains to leave the system 
through the open boundaries. This combination of slow driving and dissipation drives 
the system to the critical point of an underlying absorbing phase transition, taking place 
at a critical value £ c of the grain density £ = N/ L d [36l [37] . 

Here, we consider the two dimensional version of the model, and study the time 
series V(t) measuring the number of toppling events as a function of time, with one 
parallel update of the lattice defining the unit of time. Fig. [9] shows an example of a 
time series V(t) from the 2c? Manna model with L = 1024. By driving the system with a 
slow constant rate, the waiting times between avalanches without any thresholding will 
follow exponential statistics. The non-thresholded avalanche durations are distributed 
according to a power law with an exponent tt ~ 1.5. The average avalanche size (s(T)) 
scales with the avalanche duration as (s(T)) ~ T 7st , with 7 st « 1.77 [10]. As the 
threshold level Vth is increased from zero, a power law part starts to emerge to the 
internal waiting time distribution, with a cut-off f* scaling roughly as 

ft ~ - Vtf, (15) 

The exponent r w assumes a value close to the one observed for the BTW model [26J, i.e. 
t w ~ 1.62 ± 0.05, see Fig. [101 Similar power law scaling is observed by considering the 
ensemble of all waiting times, f = f e U fj, where f e denotes the inter- avalanche waiting 
times. The power law exponent remains unchanged from the case where only fj's are 
considered, but the inclusion of f e 's has the effect of producing a "bump" -like cut- 
off to the distributions. Interestingly, Fig. [TTJ shows that also the duration distribution 
exponent appears to evolve towards a similar value as the threshold level V t h is increased, 
in agreement with the observations in the BTW model |26j . 




Figure 10. Top: The distributions of the intra-avalanche waiting times t% in the d = 2 
Manna model with L = 1024 as a function of the threshold level Vth- The solid line 
corresponds to a power law with an exponent t w = 1.62. The inset displays a data 
collapse of the distributions for Vth > 40 showing that the cut-off scales as Tq ~ V t ]^ 3 . 
Bottom: The distributions of all (intra-avalanche as well as inter-avalanche) waiting 
times f. The solid line corresponds to a power law with an exponent t w = 1.61. 




Figure 11. The distributions of avalanche durations T in the d = 2 Manna model 
with L — 1024 as a function of the threshold level V t h- The exponent appears to 
evolve from 1.5 to 1.63 as the threshold Vth is increased. 
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Figure 12. The distributions of waiting times f in the d — 2 Manna model with 
L = 1024 and Vth = 640, for two different strengths a of the additive Gaussian white 
noise. The exponent t w appears to change to t w ss 2.0 for short waiting times, while for 
waiting times longer than a tr-dependent crossover scale, the noise free value t w « 1.62 
persists. 

Finally, we consider the effect of noise on the observed waiting time statistics in the 
Manna model. By adding Gaussian white noise with zero mean and standard deviation 
o to the V(t)-signal, we observe a similar change in the r^-exponent as in the case of 
the simple random walk model: it assumes a value close to r w « 2 for waiting times 
smaller than some crossover scale growing with a . For longer waiting times, the scaling 
is the same as in the absence of noise, r w m 1.65, see Fig. [121 

4. Conclusions 

The simple random walk model studied here presents a simple and transparent 
illustration of the mechanism leading to power law distributions of quiet times between 
consecutive subavalanches emerging from the thresholding process. A key observation 
here is that this can be understood to follow from the symmetry between the random 
walk excursions above and below the imposed threshold. The numerical result that 
such a symmetry appears to be valid also in the non-trivial case of the Manna sandpile 
model is an interesting observation that would deserve further attention. Also there is 
an apparent difference in the avalanche duration distribution exponents in the sandpile 
models between the unthresholded and thresholded cases, which is intriguing and not 
understood. 

The other interesting observation is that additive Gaussian white noise can have 
an effect on the observed waiting time statistics. The apparent change in the power law 
exponent from the noise-free value to a value close to -2.0 both in the random walk and 
Manna models for high enough noise strength remains to be explained. This observation 
may be relevant when noisy experimental signals have to be thresholded. 

Another technique that we have used to study the effect of thresholding on 
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these artificial signals is the Power Spectrum (PS) analysis. For the signal composed 
by random walk excursions the global behaviour of the PS does not change when 
thresholding due to the RW translational invariance, obtaining P(f) = A/f 2 . The 
constant A is decreasing when the threshold is increased due to the fact that the 
total energy of the signal decreases. Similar observations, though less conclusive, can 
be obtained for the Manna model: the exponent of the power spectrum does not 
seem to change significantly from the non-thresholded case where P(f) = A/f a , with 
a = j st ~ 1.77 [10J. The exponent j st is related to the avalanche distribution exponents 
through 7 si = {j T — l)/(r s — 1). As we observe an apparent change in the value of the 
r^-exponent, the r s -exponent should thus experience a corresponding change to keep 
•y st constant. Given the above values for r T and •y st , one would thus expect r s « 1.36, 
to be compared with r ps 1.28 for V t h = [38]. However, due to the limited accuracy 
of the numerical values of these exponents, and to the small (if any) change to the PS 
exponent, we do not pursue this issue. 

To discuss the implications of the toy models we have explored we note that power- 
law waiting times (and correlations) arise only after thresholding (single) avalanches. 
Then the self-affine character of the single avalanches produces power-law (PL) waiting 
times with a particular, model-dependent exponent. For empirical examples such as 
earthquakes the question is how to interpret the PL statistics and correlations in the 
cascades of aftershocks. Based on toy models one can envision three reasons for PL 
waiting times, or two additional ones beyond such thresholding. First, the system 
would be driven externally in a non-Poissonian way such that the PL statistics ensues. 
In the case of solar physics this has been argued to be the case for the magnetosphere, 
since it might be so that the solar wind which drives it has such correlations [39|. 
Second, the models we use could be far from the real dynamics, which may be much 
more complicated. We represent (in the Manna case, for instance) the dynamics with 
a projection to a one-dimensional signal V(t), but in reality such a trick may produce 
a V(t) that appears to have correlations among avalanches while in fact the system 
evolves dynamically also when V(t)=0. A similar mechanism would be that the system 
has a memory which e.g. is influenced by the size of an avalanche and influences the 
rate at which subsequent avalanches are triggered. The usual (e.g. SOC) models do not 
exhibit such, but this does not logically exclude the possibility at all. 

Finally we point out that in addition to the temporal clustering of the avalanches, 
the effect of a finite detection threshold on the often observed spatial avalanche clustering 
could be considered as well. It is also intriguing as to what would happen to two-point 
correlations in spatio-temporal systems, in both time and space domains. 
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